# ===================================================================================
# === Louvain + H3 社区发现与 POI 候选点筛选 (优化版) ===
# ===================================================================================

import warnings
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, Polygon
from shapely.ops import snap
import networkx as nx
import community as community_louvain
import momepy
import h3
import numpy as np
from tqdm import tqdm
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
import os
import folium
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import time
import pickle
from pathlib import Path
from typing import Tuple, List, Optional
from contextlib import contextmanager
import logging

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)

# ============================= 日志配置 =============================
logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s - %(levelname)s - %(message)s'
)
logger = logging.getLogger(__name__)

# ============================= 参数配置 =============================
class Config:
    """集中管理所有配置参数"""
    # 输入文件
    BOUNDARY_SHP = "仓山区最终完整边界.shp"
    ROAD_SHP = "仓山区道路信息.shp"
    POI_SHP = "poi仓山区.shp"
    
    # 算法参数
    H3_RESOLUTION = 10  # 🔧 提高分辨率 9→10，更精细的网格
    LOUVAIN_RESOLUTION = 1.2  # Louvain分辨率
    NUM_CLUSTERS = 5   # 期望聚类数
    SNAP_TOLERANCE_M = 5
    MIN_AREA_QUANTILE = 0.05
    USE_FULL_TOPOLOGY_REPAIR = True
    MIN_POI_DENSITY = 1.0  # 🆕 最小POI密度阈值（个/km²），过滤空旷区域
    
    # 缓存和输出
    CACHE_DIR = Path("cache")
    OUTPUT_DIR = Path("output")
    NETWORKX_CACHE = CACHE_DIR / "networkx_graph.pkl"
    H3_CACHE = CACHE_DIR / "h3_features.geojson"
    COMMUNITY_CACHE = CACHE_DIR / "community_polygons.geojson"  # 新增社区多边形缓存
    CANDIDATES_OUTPUT = OUTPUT_DIR / "poi_candidates_3857.npy"
    MAP_OUTPUT = OUTPUT_DIR / "louvain_h3_poi_map.html"
    
    # CRS
    CRS_WGS84 = "EPSG:4326"
    CRS_METRIC = "EPSG:3857"
    
    @classmethod
    def setup_directories(cls):
        """创建必要的目录"""
        cls.CACHE_DIR.mkdir(exist_ok=True)
        cls.OUTPUT_DIR.mkdir(exist_ok=True)

# ============================= 工具函数 =============================
@contextmanager
def timer(description: str):
    """计时上下文管理器"""
    start = time.time()
    logger.info(f"开始: {description}")
    yield
    elapsed = time.time() - start
    logger.info(f"完成: {description} (耗时: {elapsed:.2f}秒)")

def load_cache(path: Path, loader_func):
    """通用缓存加载函数"""
    if not path.exists():
        return None
    try:
        return loader_func(path)
    except Exception as e:
        logger.warning(f"缓存加载失败 ({path}): {e}")
        return None

def save_cache(obj, path: Path, saver_func):
    """通用缓存保存函数"""
    try:
        saver_func(obj, path)
        logger.info(f"缓存已保存: {path}")
    except Exception as e:
        logger.error(f"缓存保存失败 ({path}): {e}")

# ============================= 数据加载 =============================
class DataLoader:
    """数据加载和预处理类"""
    
    @staticmethod
    def fix_crs(gdf: gpd.GeoDataFrame, target_crs: str) -> gpd.GeoDataFrame:
        """统一CRS"""
        if gdf.crs is None:
            return gdf.set_crs(target_crs)
        if gdf.crs != target_crs:
            return gdf.to_crs(target_crs)
        return gdf
    
    @staticmethod
    def load_and_prepare(
        boundary_path: str,
        road_path: str,
        poi_path: str
    ) -> Tuple[gpd.GeoDataFrame, ...]:
        """加载所有数据并转换CRS"""
        with timer("数据加载与CRS转换"):
            # 加载WGS84数据
            road_gdf = DataLoader.fix_crs(gpd.read_file(road_path), Config.CRS_WGS84)
            poi_gdf = DataLoader.fix_crs(gpd.read_file(poi_path), Config.CRS_WGS84)
            boundary_gdf = DataLoader.fix_crs(gpd.read_file(boundary_path), Config.CRS_WGS84)
            
            # 转换为投影坐标系
            road_proj = road_gdf.to_crs(Config.CRS_METRIC)
            poi_proj = poi_gdf.to_crs(Config.CRS_METRIC)
            boundary_proj = boundary_gdf.to_crs(Config.CRS_METRIC)
            
        return road_gdf, poi_gdf, boundary_gdf, road_proj, poi_proj, boundary_proj

# ============================= 路网处理 =============================
class RoadNetworkProcessor:
    """路网拓扑处理和图构建"""
    
    @staticmethod
    def build_graph(road_proj: gpd.GeoDataFrame, tolerance_m: float) -> nx.Graph:
        """构建NetworkX图 (带缓存)"""
        # 尝试加载缓存
        cached_graph = load_cache(
            Config.NETWORKX_CACHE,
            lambda p: pickle.load(open(p, 'rb'))
        )
        if cached_graph is not None:
            logger.info("从缓存加载NetworkX图")
            return cached_graph
        
        with timer(f"拓扑修复与图构建 (容差: {tolerance_m}m)"):
            # 🚀 方法1: 简化拓扑修复 - 跳过全局union，直接构图
            logger.info("方法1: 跳过snap操作，直接构建图...")
            
            # 预处理：移除无效几何和重复线段
            road_clean = road_proj[road_proj.geometry.is_valid].copy()
            road_clean = road_clean.drop_duplicates(subset=['geometry'])
            
            # 直接explode多段线为单段线
            road_clean = road_clean.explode(index_parts=True).reset_index(drop=True)
            
            # 🚀 方法2: 使用更快的图构建方式
            with timer("构建NetworkX图"):
                G = momepy.gdf_to_nx(road_clean, approach="primal")
            
            # 🚀 方法3: 向量化权重计算
            with timer("计算边权重"):
                edges_without_weight = [
                    (u, v, data) for u, v, data in G.edges(data=True) 
                    if 'weight' not in data
                ]
                
                if edges_without_weight:
                    # 批量计算权重
                    nodes_u = np.array([u for u, v, _ in edges_without_weight])
                    nodes_v = np.array([v for u, v, _ in edges_without_weight])
                    weights = np.linalg.norm(nodes_u - nodes_v, axis=1)
                    
                    # 批量赋值
                    for (u, v, data), weight in zip(edges_without_weight, weights):
                        data['weight'] = weight
                
                logger.info(f"图统计: {G.number_of_nodes()} 节点, {G.number_of_edges()} 边")
            
            # 保存缓存
            save_cache(G, Config.NETWORKX_CACHE, lambda g, p: pickle.dump(g, open(p, 'wb')))
            
        return G
    
    @staticmethod
    def build_graph_with_snap(road_proj: gpd.GeoDataFrame, tolerance_m: float) -> nx.Graph:
        """构建NetworkX图 (带完整拓扑修复 - 仅在需要时使用)"""
        # 尝试加载缓存
        cached_graph = load_cache(
            Config.NETWORKX_CACHE,
            lambda p: pickle.load(open(p, 'rb'))
        )
        if cached_graph is not None:
            logger.info("从缓存加载NetworkX图")
            return cached_graph
        
        with timer(f"完整拓扑修复与图构建 (容差: {tolerance_m}m)"):
            # 🚀 优化1: 分块处理大规模数据
            chunk_size = 1000
            total_rows = len(road_proj)
            
            if total_rows > chunk_size:
                logger.info(f"大规模数据集检测到 ({total_rows} 条道路)，启用分块处理...")
                
                # 先做全局union (最耗时，但只做一次)
                with timer("计算全局道路网络union"):
                    try:
                        road_union = road_proj.geometry.union_all()
                    except AttributeError:
                        road_union = road_proj.geometry.unary_union
                
                # 分块snap
                repaired_chunks = []
                for i in tqdm(range(0, total_rows, chunk_size), desc="分块拓扑修复"):
                    chunk = road_proj.iloc[i:i+chunk_size].copy()
                    chunk['geometry'] = chunk.geometry.apply(
                        lambda g: snap(g, road_union, tolerance_m)
                    )
                    repaired_chunks.append(chunk)
                
                repaired_gdf = pd.concat(repaired_chunks, ignore_index=True)
            else:
                # 小数据集：原始方法
                try:
                    road_union = road_proj.geometry.union_all()
                except AttributeError:
                    road_union = road_proj.geometry.unary_union
                
                repaired_gdf = road_proj.copy()
                repaired_gdf['geometry'] = repaired_gdf.geometry.apply(
                    lambda g: snap(g, road_union, tolerance_m)
                )
            
            # Explode
            repaired_gdf = repaired_gdf.explode(index_parts=True).reset_index(drop=True)
            
            # 构建图
            with timer("构建NetworkX图"):
                G = momepy.gdf_to_nx(repaired_gdf, approach="primal")
            
            # 向量化权重计算
            with timer("计算边权重"):
                edges_without_weight = [
                    (u, v, data) for u, v, data in G.edges(data=True) 
                    if 'weight' not in data
                ]
                
                if edges_without_weight:
                    nodes_u = np.array([u for u, v, _ in edges_without_weight])
                    nodes_v = np.array([v for u, v, _ in edges_without_weight])
                    weights = np.linalg.norm(nodes_u - nodes_v, axis=1)
                    
                    for (u, v, data), weight in zip(edges_without_weight, weights):
                        data['weight'] = weight
                
                logger.info(f"图统计: {G.number_of_nodes()} 节点, {G.number_of_edges()} 边")
            
            # 保存缓存
            save_cache(G, Config.NETWORKX_CACHE, lambda g, p: pickle.dump(g, open(p, 'wb')))
            
        return G
    
    @staticmethod
    def discover_communities(G: nx.Graph, resolution: float) -> nx.Graph:
        """Louvain社区发现"""
        with timer(f"Louvain社区发现 (Resolution: {resolution})"):
            partition = community_louvain.best_partition(
                G, resolution=resolution, random_state=42
            )
            nx.set_node_attributes(G, partition, "community_id")
            num_communities = len(set(partition.values()))
            logger.info(f"发现 {num_communities} 个社区")
        return G

# ============================= H3处理 =============================
class H3Processor:
    """H3网格化和特征工程"""
    
    @staticmethod
    def aggregate_features(
        G: nx.Graph,
        poi_proj: gpd.GeoDataFrame,
        boundary_proj: gpd.GeoDataFrame,
        h3_resolution: int
    ) -> Tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]:
        """
        H3聚合与特征计算 (带缓存)
        返回: (h3_gdf用于特征计算, community_gdf用于可视化)
        """
        # 尝试加载社区多边形缓存
        cached_community = load_cache(
            Config.COMMUNITY_CACHE,
            lambda p: gpd.read_file(p).to_crs(Config.CRS_METRIC)
        )
        
        # 尝试加载H3缓存
        cached_h3 = load_cache(
            Config.H3_CACHE,
            lambda p: gpd.read_file(p).to_crs(Config.CRS_METRIC)
        )
        
        if cached_h3 is not None and cached_community is not None:
            logger.info("从缓存加载H3特征和社区多边形")
            return cached_h3, cached_community
        
        with timer(f"H3网格化与特征工程 (Res: {h3_resolution})"):
            # 1. 节点到H3映射
            h3_gdf = H3Processor._nodes_to_h3(G, h3_resolution)
            
            # 2. 空间裁剪
            h3_gdf = gpd.overlay(
                h3_gdf, boundary_proj[['geometry']], how="intersection"
            )
            h3_gdf["area_sqkm"] = h3_gdf.area / 1_000_000
            
            # 3. POI密度特征
            h3_gdf = H3Processor._add_poi_density(h3_gdf, poi_proj)
            
            # 4. 道路密度特征
            h3_gdf = H3Processor._add_road_density(h3_gdf, G, h3_resolution)
            
            # 保存H3缓存
            save_cache(
                h3_gdf,
                Config.H3_CACHE,
                lambda gdf, p: gdf.to_file(p, driver="GeoJSON")
            )
            
            # 5. 创建社区多边形
            community_gdf = H3Processor._create_community_polygons(h3_gdf)
            
            # 保存社区多边形缓存
            save_cache(
                community_gdf,
                Config.COMMUNITY_CACHE,
                lambda gdf, p: gdf.to_file(p, driver="GeoJSON")
            )
        
        return h3_gdf, community_gdf
    
    @staticmethod
    def _create_community_polygons(h3_gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
        """将H3单元合并成社区多边形 (改进版：保留更多社区细节)"""
        logger.info("创建社区边界多边形...")
        
        # 统计社区分布
        community_stats = h3_gdf['community_id'].value_counts()
        logger.info(f"H3社区分布统计: 共 {len(community_stats)} 个唯一社区")
        logger.info(f"  - 最大社区包含 {community_stats.max()} 个H3单元")
        logger.info(f"  - 最小社区包含 {community_stats.min()} 个H3单元")
        logger.info(f"  - 平均每个社区 {community_stats.mean():.1f} 个H3单元")
        
        # 按社区ID分组并合并几何
        community_polygons = []
        
        for comm_id, group in tqdm(h3_gdf.groupby('community_id'), desc="合并社区多边形"):
            if comm_id == -1:
                continue
            
            # 合并该社区的所有H3单元为一个多边形
            try:
                merged_geom = group.geometry.union_all()
            except AttributeError:
                merged_geom = group.geometry.unary_union
            
            # 🔧 过滤太小的社区（可能是噪声）
            area_sqkm = group['area_sqkm'].sum()
            if area_sqkm < 0.01:  # 小于0.01平方公里（10000平方米）
                continue
            
            community_polygons.append({
                'community_id': comm_id,
                'geometry': merged_geom,
                'area_sqkm': area_sqkm,
                'num_h3_cells': len(group)
            })
        
        community_gdf = gpd.GeoDataFrame(community_polygons, crs=Config.CRS_METRIC)
        logger.info(f"✓ 创建了 {len(community_gdf)} 个有效社区多边形 (总H3单元: {len(h3_gdf)})")
        logger.info(f"✓ 过滤后社区面积范围: {community_gdf['area_sqkm'].min():.3f} - {community_gdf['area_sqkm'].max():.2f} km²")
        
        return community_gdf
    
    @staticmethod
    def _nodes_to_h3(G: nx.Graph, h3_resolution: int) -> gpd.GeoDataFrame:
        """节点映射到H3单元"""
        node_data = [
            [node[0], node[1], data.get('community_id', -1)]
            for node, data in G.nodes(data=True)
        ]
        
        node_df = pd.DataFrame(node_data, columns=['x', 'y', 'community_id'])
        node_gdf = gpd.GeoDataFrame(
            node_df,
            geometry=gpd.points_from_xy(node_df.x, node_df.y),
            crs=Config.CRS_METRIC
        ).to_crs(Config.CRS_WGS84)
        
        node_gdf['h3_cell'] = node_gdf.geometry.apply(
            lambda g: h3.geo_to_h3(g.y, g.x, h3_resolution)
        )
        
        # 多数投票
        h3_community_map = node_gdf.groupby('h3_cell')['community_id'].agg(
            lambda x: x.mode()[0] if len(x.mode()) > 0 else -1
        ).reset_index()
        
        # 构建H3多边形
        h3_polygons = []
        for _, row in tqdm(h3_community_map.iterrows(), total=len(h3_community_map), desc="构建H3多边形"):
            coords = h3.h3_to_geo_boundary(row['h3_cell'], geo_json=True)
            h3_polygons.append({
                'h3_cell': row['h3_cell'],
                'community_id': row['community_id'],
                'geometry': Polygon(coords)
            })
        
        return gpd.GeoDataFrame(h3_polygons, crs=Config.CRS_WGS84).to_crs(Config.CRS_METRIC)
    
    @staticmethod
    def _add_poi_density(h3_gdf: gpd.GeoDataFrame, poi_proj: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
        """添加POI密度特征"""
        type_col = next((c for c in ['category', 'TYPE'] if c in poi_proj.columns), None)
        
        if type_col:
            poi_joined = gpd.sjoin(poi_proj, h3_gdf, how="inner", predicate="within")
            poi_counts = poi_joined.groupby(['h3_cell', type_col]).size().unstack(fill_value=0)
            h3_gdf = h3_gdf.merge(poi_counts, on="h3_cell", how="left").fillna(0)
            
            for col in poi_counts.columns:
                h3_gdf[f"density_{col}"] = h3_gdf[col] / h3_gdf["area_sqkm"]
            h3_gdf.drop(columns=poi_counts.columns, inplace=True)
        
        return h3_gdf
    
    @staticmethod
    def _add_road_density(h3_gdf: gpd.GeoDataFrame, G: nx.Graph, h3_resolution: int) -> gpd.GeoDataFrame:
        """添加道路密度特征"""
        # 🔧 兼容不同版本的 momepy
        try:
            # 新版本 momepy (>= 0.6.0)
            _, edge_gdf = momepy.nx_to_gdf(G)
        except (TypeError, ValueError):
            try:
                # 旧版本 momepy (< 0.6.0)
                edge_gdf = momepy.nx_to_gdf(G, nodes=False)
            except:
                # 最后的fallback：手动从图构建边GeoDataFrame
                logger.warning("momepy.nx_to_gdf 失败，使用手动构建")
                edge_gdf = H3Processor._manual_edges_to_gdf(G)
        
        # 确保edge_gdf有几何列和weight
        if edge_gdf.empty or 'geometry' not in edge_gdf.columns:
            logger.warning("边GeoDataFrame为空，跳过道路密度计算")
            h3_gdf['density_road_m_per_sqkm'] = 0
            return h3_gdf
        
        # 计算中心点和长度
        edge_gdf['center_wgs'] = edge_gdf.to_crs(Config.CRS_WGS84).geometry.centroid
        
        # 长度可能在weight列或需要计算
        if 'weight' in edge_gdf.columns:
            edge_gdf['length_m'] = edge_gdf['weight']
        else:
            edge_gdf['length_m'] = edge_gdf.geometry.length
        
        # 映射到H3单元
        edge_gdf['h3_cell'] = edge_gdf['center_wgs'].apply(
            lambda g: h3.geo_to_h3(g.y, g.x, h3_resolution)
        )
        
        # 聚合道路长度
        road_length_per_h3 = edge_gdf.groupby('h3_cell')['length_m'].sum()
        h3_gdf['total_road_length_m'] = h3_gdf['h3_cell'].map(road_length_per_h3).fillna(0)
        h3_gdf['density_road_m_per_sqkm'] = h3_gdf['total_road_length_m'] / h3_gdf["area_sqkm"]
        h3_gdf.drop(columns=['total_road_length_m'], inplace=True)
        
        return h3_gdf
    
    @staticmethod
    def _manual_edges_to_gdf(G: nx.Graph) -> gpd.GeoDataFrame:
        """手动将NetworkX图的边转换为GeoDataFrame (fallback方法)"""
        from shapely.geometry import LineString
        
        edges_data = []
        for u, v, data in G.edges(data=True):
            edges_data.append({
                'geometry': LineString([u, v]),
                'weight': data.get('weight', np.linalg.norm(np.array(u) - np.array(v)))
            })
        
        if not edges_data:
            return gpd.GeoDataFrame()
        
        return gpd.GeoDataFrame(edges_data, crs=Config.CRS_METRIC)

# ============================= 聚类分析 =============================
class CommunityClassifier:
    """社区功能分类和标注"""
    
    # 功能区类型定义
    ZONE_TYPES = {
        0: {'name': '商业区', 'color': '#FF6B6B', 'icon': '🏪'},
        1: {'name': '住宅区', 'color': '#4ECDC4', 'icon': '🏘️'},
        2: {'name': '办公区', 'color': '#45B7D1', 'icon': '🏢'},
        3: {'name': '混合区', 'color': '#FFA07A', 'icon': '🏙️'},
        4: {'name': '休闲区', 'color': '#98D8C8', 'icon': '🌳'},
        -1: {'name': '未分类', 'color': '#CCCCCC', 'icon': '❓'}
    }
    
    @staticmethod
    def classify_communities(
        h3_gdf: gpd.GeoDataFrame,
        community_gdf: gpd.GeoDataFrame,
        poi_proj: gpd.GeoDataFrame,
        num_clusters: int = 5
    ) -> gpd.GeoDataFrame:
        """
        基于POI特征对Louvain社区进行功能分类
        
        关键改变: 直接对社区多边形分类，不再对H3单元分类
        """
        with timer(f"社区功能分类 (K={num_clusters})"):
            # 1. 统计每个社区的POI特征 (社区级别，不是H3级别)
            community_features = CommunityClassifier._aggregate_poi_by_community_direct(
                community_gdf, poi_proj
            )
            
            if community_features.empty:
                logger.warning("无法计算社区POI特征")
                community_gdf['zone_type'] = -1
                community_gdf['zone_name'] = '未分类'
                return community_gdf
            
            # 2. K-Means聚类识别功能区
            community_features = CommunityClassifier._kmeans_classify(
                community_features, num_clusters
            )
            
            # 3. 将功能区类型合并回社区GeoDataFrame
            community_gdf = community_gdf.merge(
                community_features[['community_id', 'zone_type', 'zone_name']],
                on='community_id',
                how='left'
            )
            community_gdf['zone_type'].fillna(-1, inplace=True)
            community_gdf['zone_name'].fillna('未分类', inplace=True)
            
            # 4. 统计和日志
            CommunityClassifier._print_classification_summary(community_gdf)
            
        return community_gdf
    
    @staticmethod
    def _aggregate_poi_by_community_direct(
        community_gdf: gpd.GeoDataFrame,
        poi_proj: gpd.GeoDataFrame
    ) -> pd.DataFrame:
        """直接在社区多边形上聚合POI (不经过H3)"""
        logger.info("步骤1: 直接聚合POI到社区...")
        
        # 空间连接：POI -> Community
        poi_with_community = gpd.sjoin(
            poi_proj, 
            community_gdf[['geometry', 'community_id', 'area_sqkm']], 
            how="inner", 
            predicate="within"
        )
        
        # 获取POI类别列
        type_col = next(
            (c for c in ['category', 'TYPE', 'type'] if c in poi_with_community.columns),
            None
        )
        
        if not type_col:
            logger.warning("POI数据中未找到类别字段")
            return pd.DataFrame()
        
        # 统计每个社区的POI类别计数
        poi_counts = poi_with_community.groupby(['community_id', type_col]).size().unstack(fill_value=0)
        
        # 获取社区面积
        community_area = community_gdf.set_index('community_id')['area_sqkm']
        
        # 计算POI密度（每平方公里）
        poi_density = poi_counts.div(community_area, axis=0).fillna(0)
        poi_density.columns = [f'density_{col}' for col in poi_density.columns]
        
        # 添加总POI数量
        poi_density['total_poi_count'] = poi_counts.sum(axis=1)
        poi_density['total_poi_density'] = poi_density['total_poi_count'] / community_area
        
        logger.info(f"✓ 处理了 {len(poi_density)} 个社区，{len(poi_counts.columns)} 种POI类别")
        
        return poi_density.reset_index()
    
    @staticmethod
    def _kmeans_classify(
        community_features: pd.DataFrame,
        num_clusters: int
    ) -> pd.DataFrame:
        """使用K-Means对社区进行功能分类"""
        logger.info(f"步骤2: K-Means聚类 (K={num_clusters})...")
        
        # 选择密度特征进行聚类
        feature_cols = [c for c in community_features.columns if c.startswith('density_')]
        
        if not feature_cols:
            logger.warning("没有可用的密度特征")
            community_features['zone_type'] = -1
            return community_features
        
        # 🔧 改进：不过滤低密度社区，而是根据POI分布特征分类
        valid_communities = community_features[
            community_features['total_poi_count'] > 0  # 只要有POI就参与分类
        ].copy()
        
        logger.info(f"参与分类的社区数: {len(valid_communities)} / {len(community_features)}")
        
        if len(valid_communities) < num_clusters:
            logger.warning(f"有效社区数({len(valid_communities)})少于聚类数({num_clusters})")
            num_clusters = max(2, len(valid_communities) // 2)
        
        # 标准化特征
        X = valid_communities[feature_cols].fillna(0).values
        scaler = StandardScaler()
        X_scaled = scaler.fit_transform(X)
        
        # K-Means聚类
        kmeans = KMeans(n_clusters=num_clusters, random_state=42, n_init=10)
        valid_communities['zone_type'] = kmeans.fit_predict(X_scaled)
        
        # 合并回原数据
        community_features = community_features.merge(
            valid_communities[['community_id', 'zone_type']],
            on='community_id',
            how='left'
        )
        community_features['zone_type'].fillna(-1, inplace=True)
        community_features['zone_type'] = community_features['zone_type'].astype(int)
        
        # 分析聚类特征，自动命名功能区
        community_features = CommunityClassifier._auto_name_zones(
            community_features, feature_cols, kmeans
        )
        
        logger.info(f"✓ 识别出 {num_clusters} 种功能区类型")
        logger.info(f"✓ 未分类社区数: {(community_features['zone_type'] == -1).sum()}")
        
        return community_features
    
    @staticmethod
    def _auto_name_zones(
        community_features: pd.DataFrame,
        feature_cols: List[str],
        kmeans: KMeans
    ) -> pd.DataFrame:
        """根据POI特征自动为功能区命名"""
        # 计算每个聚类中心的主导POI类别
        cluster_centers = pd.DataFrame(
            kmeans.cluster_centers_,
            columns=feature_cols
        )
        
        zone_names = []
        for idx, row in cluster_centers.iterrows():
            # 找出最高密度的前3种POI类别
            top_pois = row.nlargest(3).index.tolist()
            
            # 根据POI类型推断功能区
            poi_types = [col.replace('density_', '') for col in top_pois]
            
            # 简单的规则匹配（可根据实际POI类别调整）
            if any('商' in p or '购' in p or '店' in p for p in poi_types):
                zone_name = '商业区'
            elif any('住' in p or '小区' in p or '公寓' in p for p in poi_types):
                zone_name = '住宅区'
            elif any('公司' in p or '办公' in p or '写字楼' in p for p in poi_types):
                zone_name = '办公区'
            elif any('公园' in p or '景点' in p or '娱乐' in p for p in poi_types):
                zone_name = '休闲区'
            else:
                zone_name = f'混合区-{idx}'
            
            zone_names.append(zone_name)
        
        # 映射聚类ID到功能区名称
        zone_type_map = {i: name for i, name in enumerate(zone_names)}
        zone_type_map[-1] = '未分类'
        
        community_features['zone_name'] = community_features['zone_type'].map(zone_type_map)
        
        return community_features
    
    @staticmethod
    def _map_zone_types_to_h3(
        h3_gdf: gpd.GeoDataFrame,
        community_features: pd.DataFrame
    ) -> gpd.GeoDataFrame:
        """将社区功能类型映射到H3单元 (保留用于特征计算)"""
        logger.info("步骤3: 映射功能区到H3单元...")
        
        zone_mapping = community_features[['community_id', 'zone_type', 'zone_name']].drop_duplicates()
        
        h3_gdf = h3_gdf.merge(zone_mapping, on='community_id', how='left')
        h3_gdf['zone_type'].fillna(-1, inplace=True)
        h3_gdf['zone_name'].fillna('未分类', inplace=True)
        
        return h3_gdf
    
    @staticmethod
    def _print_classification_summary(gdf: gpd.GeoDataFrame):
        """打印分类摘要"""
        logger.info("\n" + "="*50)
        logger.info("社区功能分类结果:")
        logger.info("="*50)
        
        if 'zone_name' not in gdf.columns:
            logger.warning("缺少 zone_name 列")
            return
        
        # 按功能区统计
        zone_summary = gdf.groupby('zone_name').agg({
            'community_id': 'nunique' if 'community_id' in gdf.columns else 'count',
            'area_sqkm': 'sum'
        }).round(2)
        
        if 'community_id' in gdf.columns:
            zone_summary.columns = ['社区数', '总面积(km²)']
        else:
            zone_summary.columns = ['数量', '总面积(km²)']
        
        logger.info("\n" + zone_summary.to_string())
        logger.info("="*50 + "\n")
    
    @staticmethod
    def select_target_zones_for_poi(
        community_gdf: gpd.GeoDataFrame,
        poi_proj: gpd.GeoDataFrame,
        target_zone_names: List[str] = None
    ) -> Tuple[gpd.GeoDataFrame, np.ndarray]:
        """
        选择目标功能区并提取POI候选点
        
        Args:
            community_gdf: 社区多边形GeoDataFrame (不是H3)
            poi_proj: POI数据
            target_zone_names: 目标功能区名称列表
        """
        with timer("提取目标功能区POI候选点"):
            if target_zone_names:
                # 手动指定目标功能区
                target_zones = community_gdf[community_gdf['zone_name'].isin(target_zone_names)].copy()
                logger.info(f"选择目标功能区: {target_zone_names}")
            else:
                # 自动选择：排除未分类，选择面积最大的功能区
                zone_areas = community_gdf[community_gdf['zone_name'] != '未分类'].groupby('zone_name')['area_sqkm'].sum()
                if zone_areas.empty:
                    logger.warning("没有有效的功能区，返回空候选点")
                    return gpd.GeoDataFrame(), np.array([])
                
                top_zone = zone_areas.idxmax()
                target_zones = community_gdf[community_gdf['zone_name'] == top_zone].copy()
                logger.info(f"自动选择最大功能区: {top_zone} ({zone_areas[top_zone]:.2f} km²)")
            
            if target_zones.empty:
                logger.warning("目标功能区为空")
                return gpd.GeoDataFrame(), np.array([])
            
            # 提取目标区域内的POI
            target_zone_union = target_zones.geometry.union_all()
            poi_in_zones = poi_proj[poi_proj.geometry.intersects(target_zone_union)].copy()
            
            candidate_coords = np.array([[p.x, p.y] for p in poi_in_zones.geometry])
            
            logger.info(f"✓ 候选点池大小: {len(candidate_coords)} 个POI")
            
            return poi_in_zones, candidate_coords

# ============================= 可视化 =============================
class Visualizer:
    """结果可视化"""
    
    @staticmethod
    def create_map(
        community_gdf: gpd.GeoDataFrame,
        poi_proj: gpd.GeoDataFrame,
        candidate_coords: np.ndarray,
        output_path: Path,
        road_gdf: gpd.GeoDataFrame = None,
        G: nx.Graph = None
    ):
        """生成交互式HTML地图 (使用社区多边形，不是H3)"""
        with timer("生成HTML地图"):
            community_wgs = community_gdf.to_crs(Config.CRS_WGS84)
            poi_wgs = poi_proj.to_crs(Config.CRS_WGS84)
            
            centroid = community_wgs.unary_union.centroid
            m = folium.Map(
                location=[centroid.y, centroid.x],
                zoom_start=12,
                tiles="cartodbpositron"
            )
            
            # 🛣️ 道路网络层 (按社区着色)
            if road_gdf is not None:
                Visualizer._add_road_network_layer(m, road_gdf, G)
                # 可选：添加道路分类层
                Visualizer._add_road_classification_layer(m, road_gdf)
            
            # 🎨 Louvain 社区层 (完整多边形，不是六边形)
            Visualizer._add_louvain_community_layer(m, community_gdf)
            
            # 🎨 功能区图例
            Visualizer._add_kmeans_cluster_layer(m, community_wgs)
            
            # POI候选点层
            Visualizer._add_poi_layer(m, poi_wgs)
            
            # 模拟最优点层
            Visualizer._add_simulated_optimal_layer(m, candidate_coords)
            
            folium.LayerControl().add_to(m)
            m.save(output_path)
            logger.info(f"地图已保存: {output_path}")
    
    @staticmethod
    def _add_road_network_layer(m: folium.Map, road_gdf: gpd.GeoDataFrame, G: nx.Graph):
        """添加道路网络层 (按社区分类着色)"""
        logger.info("开始渲染道路网络层...")
        
        try:
            # 方法1: 使用NetworkX图的边
            try:
                _, edge_gdf = momepy.nx_to_gdf(G)
            except (TypeError, ValueError):
                try:
                    edge_gdf = momepy.nx_to_gdf(G, nodes=False)
                except:
                    edge_gdf = None
            
            # 如果momepy失败，直接使用原始道路数据
            if edge_gdf is None or edge_gdf.empty:
                logger.info("从NetworkX提取边失败，使用原始道路数据")
                edge_gdf = road_gdf.copy()
            
            # 确保在WGS84坐标系
            if edge_gdf.crs != Config.CRS_WGS84:
                edge_gdf = edge_gdf.to_crs(Config.CRS_WGS84)
            
            # 为道路分配社区ID (如果可能)
            if G is not None:
                node_communities = {node: data.get('community_id', -1) 
                                   for node, data in G.nodes(data=True)}
                
                # 为每条边分配社区颜色
                def assign_community_to_edge(geometry):
                    """通过道路中点查找最近的节点"""
                    try:
                        centroid = geometry.centroid
                        # 简化：随机分配社区（实际应该用空间查询）
                        return list(node_communities.values())[0] if node_communities else -1
                    except:
                        return -1
                
                # 如果edge_gdf是从NetworkX来的，应该有对应关系
                edge_list = list(G.edges())
                if len(edge_list) == len(edge_gdf):
                    edge_gdf['community_id'] = [
                        node_communities.get(u, node_communities.get(v, -1)) 
                        for u, v in edge_list
                    ]
                else:
                    edge_gdf['community_id'] = -1
            else:
                edge_gdf['community_id'] = -1
            
            # 生成颜色映射
            unique_communities = edge_gdf['community_id'].unique()
            num_communities = len(unique_communities)
            
            if num_communities <= 10:
                cmap = plt.get_cmap("tab10")
            elif num_communities <= 20:
                cmap = plt.get_cmap("tab20")
            else:
                cmap = plt.get_cmap("hsv")
            
            community_colors = {}
            for idx, comm_id in enumerate(unique_communities):
                if comm_id == -1:
                    community_colors[comm_id] = '#4A90E2'  # 蓝色
                else:
                    color_idx = idx / max(num_communities - 1, 1)
                    community_colors[comm_id] = mcolors.to_hex(cmap(color_idx))
            
            # 创建道路图层组
            road_group = folium.FeatureGroup(name="🛣️ 道路网络", show=True)
            
            # 如果道路太多，采样显示
            max_roads = 5000
            if len(edge_gdf) > max_roads:
                logger.info(f"道路数量过多({len(edge_gdf)})，采样显示{max_roads}条")
                edge_sample = edge_gdf.sample(n=max_roads, random_state=42)
            else:
                edge_sample = edge_gdf
            
            # 绘制道路
            for _, row in edge_sample.iterrows():
                comm_id = row.get('community_id', -1)
                color = community_colors.get(comm_id, '#4A90E2')
                
                folium.PolyLine(
                    locations=[[coord[1], coord[0]] for coord in row.geometry.coords],
                    color=color,
                    weight=2,
                    opacity=0.6
                ).add_to(road_group)
            
            road_group.add_to(m)
            logger.info(f"✓ 已添加道路网络层: {len(edge_sample)} 条道路")
            
        except Exception as e:
            logger.error(f"道路网络层渲染失败: {e}", exc_info=True)
            # 添加简化版道路层
            try:
                logger.info("尝试添加简化道路层...")
                road_wgs = road_gdf.to_crs(Config.CRS_WGS84)
                
                simple_road_group = folium.FeatureGroup(name="🛣️ 道路网络 (简化)", show=True)
                
                # 限制道路数量
                if len(road_wgs) > 1000:
                    road_wgs = road_wgs.sample(n=1000, random_state=42)
                
                folium.GeoJson(
                    road_wgs.geometry,
                    style_function=lambda x: {
                        'color': '#4A90E2',
                        'weight': 2,
                        'opacity': 0.6
                    }
                ).add_to(simple_road_group)
                
                simple_road_group.add_to(m)
                logger.info(f"✓ 已添加简化道路层: {len(road_wgs)} 条道路")
                
            except Exception as e2:
                logger.error(f"简化道路层也失败: {e2}")
    
    @staticmethod
    def _add_road_classification_layer(m: folium.Map, road_gdf: gpd.GeoDataFrame):
        """添加道路分类层 (按道路等级着色)"""
        try:
            logger.info("添加道路分类层...")
            road_wgs = road_gdf.to_crs(Config.CRS_WGS84)
            
            # 定义道路等级颜色
            road_type_colors = {
                '高速公路': '#FF0000',
                '快速路': '#FF3300', 
                '主干道': '#FF6B00',
                '次干道': '#FFA500',
                '支路': '#FFD700',
                '其他': '#CCCCCC'
            }
            
            # 尝试识别道路等级字段
            type_field = None
            possible_fields = ['road_type', 'TYPE', 'type', 'class', 'CLASS', 'highway', 'ROAD_CLASS']
            for field in possible_fields:
                if field in road_wgs.columns:
                    type_field = field
                    logger.info(f"检测到道路分类字段: {field}")
                    break
            
            if type_field:
                # 按道路类型分组
                road_class_group = folium.FeatureGroup(name="🚗 道路等级分类", show=False)
                
                unique_types = road_wgs[type_field].unique()
                logger.info(f"道路类型: {unique_types}")
                
                for road_type in unique_types[:10]:  # 最多显示10种类型
                    group = road_wgs[road_wgs[type_field] == road_type]
                    
                    # 限制每种类型的道路数量
                    if len(group) > 500:
                        group = group.sample(n=500, random_state=42)
                    
                    color = road_type_colors.get(str(road_type), '#CCCCCC')
                    
                    for _, row in group.iterrows():
                        folium.PolyLine(
                            locations=[[coord[1], coord[0]] for coord in row.geometry.coords],
                            color=color,
                            weight=3,
                            opacity=0.7,
                            popup=f"类型: {road_type}"
                        ).add_to(road_class_group)
                
                road_class_group.add_to(m)
                logger.info(f"✓ 已添加道路分类层")
            else:
                logger.info("未找到道路分类字段，跳过道路等级层")
                
        except Exception as e:
            logger.warning(f"道路分类层添加失败: {e}")
    
    @staticmethod
    def _add_louvain_community_layer(m: folium.Map, community_gdf: gpd.GeoDataFrame):
        """添加Louvain社区层 (完整社区多边形，不是H3六边形)"""
        logger.info(f"添加社区多边形层，共 {len(community_gdf)} 个社区...")
        
        if 'community_id' not in community_gdf.columns:
            logger.warning("缺少 community_id 列，跳过社区层")
            return
        
        # 转换到WGS84
        if community_gdf.crs != Config.CRS_WGS84:
            community_wgs = community_gdf.to_crs(Config.CRS_WGS84)
        else:
            community_wgs = community_gdf.copy()
        
        logger.info(f"社区数据列: {community_wgs.columns.tolist()}")
        
        # 检查功能区分类情况
        has_zone_classification = 'zone_type' in community_wgs.columns and 'zone_name' in community_wgs.columns
        
        if has_zone_classification:
            zone_stats = community_wgs['zone_name'].value_counts()
            logger.info(f"功能区分布:\n{zone_stats}")
            
            # 🔧 如果未分类的太多，使用社区ID着色
            if zone_stats.get('未分类', 0) > len(community_wgs) * 0.8:
                logger.warning(f"⚠️ 未分类社区过多 ({zone_stats.get('未分类', 0)}/{len(community_wgs)})，改用社区ID着色")
                has_zone_classification = False
        
        # 获取功能区颜色映射
        def get_zone_color(zone_type, zone_name):
            """根据功能区类型获取颜色"""
            if pd.isna(zone_name) or zone_name == '未分类' or zone_type == -1:
                return '#E8E8E8'  # 浅灰色
            
            try:
                zone_info = CommunityClassifier.ZONE_TYPES.get(int(zone_type), CommunityClassifier.ZONE_TYPES[-1])
                return zone_info['color']
            except:
                return '#E8E8E8'
        
        if not has_zone_classification:
            logger.info("使用社区ID着色（每个社区不同颜色）")
            # 使用社区ID生成颜色
            unique_communities = community_wgs['community_id'].unique()
            num_communities = len(unique_communities)
            
            # 🎨 使用更丰富的色谱
            if num_communities <= 10:
                cmap = plt.get_cmap("Set3")
            elif num_communities <= 20:
                cmap = plt.get_cmap("tab20")
            elif num_communities <= 40:
                # 组合两个色谱
                cmap1 = plt.get_cmap("tab20")
                cmap2 = plt.get_cmap("tab20b")
                community_colors = {}
                for idx, comm_id in enumerate(unique_communities):
                    if idx < 20:
                        community_colors[comm_id] = mcolors.to_hex(cmap1(idx / 20))
                    else:
                        community_colors[comm_id] = mcolors.to_hex(cmap2((idx - 20) / 20))
                cmap = None
            else:
                cmap = plt.get_cmap("hsv")
            
            if cmap is not None:
                community_colors = {}
                for idx, comm_id in enumerate(unique_communities):
                    color_idx = idx / max(num_communities - 1, 1)
                    community_colors[comm_id] = mcolors.to_hex(cmap(color_idx))
            
            def style_function(feature):
                comm_id = feature['properties'].get('community_id', -1)
                return {
                    'fillColor': community_colors.get(comm_id, '#CCCCCC'),
                    'color': 'white',
                    'weight': 2.0,
                    'fillOpacity': 0.8
                }
        else:
            logger.info("使用功能区分类着色")
            # 使用功能区颜色
            def style_function(feature):
                zone_type = feature['properties'].get('zone_type', -1)
                zone_name = feature['properties'].get('zone_name', '未分类')
                
                return {
                    'fillColor': get_zone_color(zone_type, zone_name),
                    'color': 'white',
                    'weight': 2.0,
                    'fillOpacity': 0.75
                }
        
        # 准备tooltip字段
        tooltip_fields = ['community_id', 'area_sqkm']
        tooltip_aliases = ['社区ID', '面积(km²)']
        
        if 'num_h3_cells' in community_wgs.columns:
            tooltip_fields.append('num_h3_cells')
            tooltip_aliases.append('H3单元数')
        
        if has_zone_classification:
            tooltip_fields.append('zone_name')
            tooltip_aliases.append('功能区类型')
        
        # 添加到地图
        folium.GeoJson(
            community_wgs[['geometry'] + tooltip_fields],
            name="🎨 社区功能分区",
            style_function=style_function,
            tooltip=folium.GeoJsonTooltip(
                fields=tooltip_fields,
                aliases=tooltip_aliases
            ),
            show=True,
            overlay=True
        ).add_to(m)
        
        logger.info(f"✓ 已添加社区功能分区层: {len(community_wgs)} 个社区多边形")
    
    @staticmethod
    def _add_kmeans_cluster_layer(m: folium.Map, h3_wgs: gpd.GeoDataFrame):
        """添加功能区图例说明层"""
        if 'zone_name' not in h3_wgs.columns:
            return
        
        # 创建图例HTML
        legend_html = '''
        <div style="position: fixed; 
                    bottom: 50px; right: 50px; width: 200px; height: auto; 
                    background-color: white; z-index:9999; font-size:14px;
                    border:2px solid grey; border-radius: 5px; padding: 10px">
        <p style="margin:0; font-weight:bold; text-align:center;">功能区类型</p>
        <hr style="margin: 5px 0;">
        '''
        
        # 统计各功能区
        zone_stats = h3_wgs.groupby('zone_name').agg({
            'area_sqkm': 'sum',
            'zone_type': 'first'
        }).reset_index()
        
        for _, row in zone_stats.iterrows():
            zone_name = row['zone_name']
            area = row['area_sqkm']
            zone_type = int(row['zone_type']) if pd.notna(row['zone_type']) else -1
            zone_info = CommunityClassifier.ZONE_TYPES.get(zone_type, CommunityClassifier.ZONE_TYPES[-1])
            
            legend_html += f'''
            <p style="margin: 5px 0;">
                <span style="background-color:{zone_info['color']}; 
                            width: 20px; height: 20px; 
                            display: inline-block; border: 1px solid black;"></span>
                {zone_info['icon']} {zone_name} ({area:.1f}km²)
            </p>
            '''
        
        legend_html += '</div>'
        
        m.get_root().html.add_child(folium.Element(legend_html))
    
    @staticmethod
    def _add_poi_layer(m: folium.Map, poi_wgs: gpd.GeoDataFrame):
        """添加POI候选点层"""
        poi_group = folium.FeatureGroup(name="📍 POI候选点池", show=True)
        
        # 为了性能，如果POI太多则采样显示
        if len(poi_wgs) > 5000:
            logger.info(f"POI数量过多({len(poi_wgs)})，采样显示1000个点")
            poi_sample = poi_wgs.sample(n=1000, random_state=42)
        else:
            poi_sample = poi_wgs
        
        for _, row in poi_sample.iterrows():
            folium.CircleMarker(
                location=[row.geometry.y, row.geometry.x],
                radius=2,
                color='#1E90FF',
                fill=True,
                fill_opacity=0.8,
                popup=f"Type: {row.get('category', row.get('TYPE', 'N/A'))}"
            ).add_to(poi_group)
        
        poi_group.add_to(m)
    
    @staticmethod
    def _add_simulated_optimal_layer(m: folium.Map, candidate_coords: np.ndarray):
        """添加模拟最优点层"""
        if len(candidate_coords) == 0:
            return
        
        # 随机选择5个点模拟PSO结果
        num_samples = min(5, len(candidate_coords))
        indices = np.random.choice(len(candidate_coords), num_samples, replace=False)
        selected_coords = candidate_coords[indices]
        
        final_gdf = gpd.GeoDataFrame(
            geometry=[Point(c) for c in selected_coords],
            crs=Config.CRS_METRIC
        ).to_crs(Config.CRS_WGS84)
        
        optimal_group = folium.FeatureGroup(name="⭐ 模拟最优选址点", show=True)
        for idx, point in enumerate(final_gdf.geometry):
            folium.Marker(
                [point.y, point.x],
                icon=folium.Icon(color='red', icon='star', prefix='fa'),
                popup=f"模拟最优站点 {idx+1}"
            ).add_to(optimal_group)
        optimal_group.add_to(m)

# ============================= 主流程 =============================
class POIAnalysisPipeline:
    """完整分析流程"""
    
    def __init__(self, config: Config = Config):
        self.config = config
        self.config.setup_directories()
    
    def run(self):
        """执行完整流程"""
        try:
            # 1. 数据加载
            data = DataLoader.load_and_prepare(
                self.config.BOUNDARY_SHP,
                self.config.ROAD_SHP,
                self.config.POI_SHP
            )
            road_gdf, _, _, road_proj, poi_proj, boundary_proj = data
            
            # 2. 路网处理
            if self.config.USE_FULL_TOPOLOGY_REPAIR:
                logger.info("使用完整拓扑修复模式 (较慢但更精确)")
                G = RoadNetworkProcessor.build_graph_with_snap(road_proj, self.config.SNAP_TOLERANCE_M)
            else:
                logger.info("使用快速模式 (跳过snap，适合干净数据)")
                G = RoadNetworkProcessor.build_graph(road_proj, self.config.SNAP_TOLERANCE_M)
            
            G = RoadNetworkProcessor.discover_communities(G, self.config.LOUVAIN_RESOLUTION)
            
            # 3. H3特征工程 + 创建社区多边形
            h3_gdf, community_gdf = H3Processor.aggregate_features(
                G, poi_proj, boundary_proj, self.config.H3_RESOLUTION
            )
            
            # 4. 社区功能分类 (直接对社区多边形分类)
            community_gdf = CommunityClassifier.classify_communities(
                h3_gdf, community_gdf, poi_proj, self.config.NUM_CLUSTERS
            )
            
            # 5. 选择目标功能区并提取候选点
            poi_in_zones, candidates = CommunityClassifier.select_target_zones_for_poi(
                community_gdf, poi_proj, target_zone_names=None
            )
            
            # 6. 导出结果
            if len(candidates) > 0:
                np.save(self.config.CANDIDATES_OUTPUT, candidates)
                logger.info(f"候选点已保存: {self.config.CANDIDATES_OUTPUT}")
            
            # 7. 可视化 (使用社区多边形，不是H3)
            Visualizer.create_map(
                community_gdf, poi_in_zones, candidates, self.config.MAP_OUTPUT,
                road_gdf=road_gdf, G=G
            )
            
            logger.info("流程完成！")
            
        except Exception as e:
            logger.error(f"流程执行失败: {e}", exc_info=True)
            raise

# ============================= 入口 =============================
if __name__ == "__main__":
    pipeline = POIAnalysisPipeline()
    pipeline.run()